Predizendo Chuva na Australia
Introdução
Os incêndios florestais na Austrália são uma ocorrência regular que teve um papel significativo na formação da natureza do continente ao longo de milhões de anos. No entanto, os incêndios podem causar danos materiais significativos e levar à perda de vidas humanas e animais. Os incêndios florestais mataram cerca de 800 pessoas e bilhões de animais desde 1851. Estima-se que entre 2019 e 2020 os incêndios mataram pelo menos 33 pessoas e mais de 3 bilhões de animais.
Geralmente, os incêndios mais devastadores são precedidos por altas temperaturas, baixa umidade relativa e fortes ventos, que criam as condições ideais para a rápida propagação do fogo. A chuva é a maior aliada ao combate ao fogo e na mitigação do surgimento de novos de focos de incêndios. Outro problema que agrava a situação na Austrália, é a forte dependência do país de civis para conter seus incêndios, cerca de 90% dos bombeiros combatendo fogos florestais são voluntários.
Portanto, considerando os recursos limitados para combate aos incêndios florestais, ainda mais com o aumento das mudanças climáticas que estão tornando esses eventos mais frequentes e devastadores, saber se irá chover em uma região específica dadas algumas medições é importante para direcionar os esforços onde eles são mais necessários.
Dataset
O conjunto de dados utilizado contém cerca de 10 anos de observações diárias tiradas de estações meteorológicas em 49 locais na Austrália. O dataset é composto por informações como Data e Local da observação, dados meteorológicos do dia no local como Direção do vento, humidade e temperatura além de duas features categóricas: RainToday e RainTomorrow que indicam se choveu 1 mm ou mais no dia da observação e no dia seguinte, respectivamente.
Visualizando os Dados
dataset = tibble(read.csv('weatherAUS.csv'))
print(dataset)## # A tibble: 145,460 x 23
## Date Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 2008-12-01 Albury 13.4 22.9 0.6 NA NA W
## 2 2008-12-02 Albury 7.4 25.1 0 NA NA WNW
## 3 2008-12-03 Albury 12.9 25.7 0 NA NA WSW
## 4 2008-12-04 Albury 9.2 28 0 NA NA NE
## 5 2008-12-05 Albury 17.5 32.3 1 NA NA W
## 6 2008-12-06 Albury 14.6 29.7 0.2 NA NA WNW
## 7 2008-12-07 Albury 14.3 25 0 NA NA W
## 8 2008-12-08 Albury 7.7 26.7 0 NA NA W
## 9 2008-12-09 Albury 9.7 31.9 0 NA NA NNW
## 10 2008-12-10 Albury 13.1 30.1 1.4 NA NA W
## # ... with 145,450 more rows, and 15 more variables: WindGustSpeed <int>,
## # WindDir9am <chr>, WindDir3pm <chr>, WindSpeed9am <int>, WindSpeed3pm <int>,
## # Humidity9am <int>, Humidity3pm <int>, Pressure9am <dbl>, Pressure3pm <dbl>,
## # Cloud9am <int>, Cloud3pm <int>, Temp9am <dbl>, Temp3pm <dbl>,
## # RainToday <chr>, RainTomorrow <chr>Usando a API do Google Maps para salvar as coordenadas geográficas de cada local de medição, como existem vários lugares com o mesmo nome em diferentes países anglófonos, é necessário concatenar o nome do país no final do local para garantir que as coordenadas estarão coerentes.
LocCount <- dataset %>% count(Location)
LocCount$Category <- paste(LocCount$Location, ", Australia")# Adding Australia to the end of location name since there are many cities with the same names in other countries
cities_df <- as.data.frame(LocCount)
locations_df <- mutate_geocode(cities_df, Category) # Getting the coordinates with the google maps API
locations <- as_tibble(locations_df)
print(locations)## # A tibble: 49 x 5
## Location n Category lon lat
## <chr> <int> <chr> <dbl> <dbl>
## 1 Adelaide 3193 Adelaide , Australia 139. -34.9
## 2 Albany 3040 Albany , Australia 118. -35.0
## 3 Albury 3040 Albury , Australia 147. -36.1
## 4 AliceSprings 3040 AliceSprings , Australia 134. -23.7
## 5 BadgerysCreek 3009 BadgerysCreek , Australia 151. -33.9
## 6 Ballarat 3040 Ballarat , Australia 144. -37.6
## 7 Bendigo 3040 Bendigo , Australia 144. -36.8
## 8 Brisbane 3193 Brisbane , Australia 153. -27.5
## 9 Cairns 3040 Cairns , Australia 146. -16.9
## 10 Canberra 3436 Canberra , Australia 149. -35.3
## # ... with 39 more rowsdf <- merge(x = dataset, y = locations, by = "Location", all = TRUE) # Outer join to combine the orignial df with the location coordinates# Creating Map Canvas
aus_coord <- geocode("Australia")
map <- get_googlemap(center = c(aus_coord$lon, aus_coord$lat), zoom = 4)Mapa com as Médias de Temperaturas Máximas de Cada Local
MTemp <- df[!is.na(df$MaxTemp), ] %>% # Removing missing data
select(Location, MaxTemp, lon, lat)
aggTemp <- aggregate(MTemp$MaxTemp, by = list(Locaton=MTemp$Location , Longitude=MTemp$lon , Latitude=MTemp$lat), FUN=mean) # Aggregating by Location and the Average MaxTemp
ggmap(map) +
geom_point(data = aggTemp,
aes(x = Longitude, y = Latitude, size = x),
color = "red", alpha = 0.3) +
theme(legend.position = "left") +
xlab("Longitude") +
ylab("Latitude") +
labs(size="Average Max. Temperature (°C)")
Mapa de Calor de Pluvialidade
RF <- df[!is.na(df$Rainfall), ] %>% # Removing missing data
select(Location, Rainfall, lon, lat)
aggRain <- aggregate(RF$Rainfall, by = list(Locaton=RF$Location , Longitude=RF$lon , Latitude=RF$lat), FUN=sum) # Aggregating by Location and the sum of Rainfall
YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")
ggmap(map) +
stat_density_2d(
data = aggRain,
aes(x = Longitude, y = Latitude, fill = stat(x)),
alpha = .1,
bins = 30,
geom = "polygon"
) +
scale_fill_gradientn(colors = YlOrBr) +
theme(legend.position = "left") +
xlab("Longitude") +
ylab("Latitude") +
labs(fill="Rainfall in mm")
Média de Precipitação por Mês
RD <- df[!is.na(df$Rainfall), ] %>% # Removing missing data
select(Date, Rainfall)
# substring(RD$Date,6,7) # Month only
RainMonth = aggregate(RD$Rainfall, by=list(Month=substring(RD$Date,6,7)), FUN=mean) # Average daily Rainfall by month
ggplot(data=RainMonth, aes(x=Month, y=x, group=1)) +
geom_line()+
geom_point() +
xlab("Month") +
ylab("Average Rainfall (mm per day)")
O gráfico acima demonstra que os meses de verão tem uma média de chuva muito maior que os meses de inverno, que tendem a ser mais secos.
Relação entre incidência solar e chuva
SunRain <- df[!is.na(df$Sunshine) & !is.na(df$Rainfall),] %>% # Removing missing data
select(Sunshine, Rainfall)
RainSun = aggregate(SunRain$Rainfall, by=list(Sunshine=SunRain$Sunshine), FUN=mean) # Average daily Rainfall by hours of sunshine
ggplot(data=RainSun, aes(x=Sunshine, y=x)) +
geom_line()+
xlab("Sunshine (hours in day)") +
ylab("Average Rainfall (mm)")
Como é possível notar, existe uma relação inversa entre a incidência solar e a quantidade de chuva.
Incidência Solar por Mês
SunMonth <- df[!is.na(df$Sunshine) & !is.na(df$Date),] %>% # Removing missing data
select(Date, Sunshine)
MonthSun = aggregate(SunMonth$Sunshine, by=list(Month=substring(SunMonth$Date,6,7)), FUN=sum) # sunlight by month
ggplot(data=MonthSun, aes(x=Month, y=x, group=1)) +
geom_line()+
geom_point() +
ylab("Sunshine (hours)") +
xlab("Month") 
Comparando o gráfico acima com o de Média de Precipitação por Mês, é reforçada a hipótese da relação inversa entre a incidência solar e a quantidade de chuva,
Precipitação em relação à direção do vento
DD<- df[!is.na(df$Rainfall), ] %>% # Removing missing data
select(Date, Rainfall,WindDir3pm, WindDir9am)
# Aggregating by Wind Direction and the avg Rainfall for 3pm and 9am
aggDir3 <- aggregate(DD$Rainfall, by = list(Dir3pm=DD$WindDir3pm), FUN=mean)
aggDir9 <- aggregate(DD$Rainfall, by = list(Dir9am=DD$WindDir9am), FUN=mean)
# Transposing the Dataframes so the columns are the mean Rainfall for each direction and the rows are the differente times
final_df3 <- as.data.frame(t(aggDir3))
names(final_df3) <- lapply(final_df3[1, ], as.character)
final_df3 <- final_df3[-1,]
rownames(final_df3) <- "3pm"
final_df9 <- as.data.frame(t(aggDir9))
names(final_df9) <- lapply(final_df9[1, ], as.character)
final_df9 <- final_df9[-1,]
rownames(final_df9) <- "9am"
final_df <- rbind(final_df3, final_df9)
final_df <- final_df %>% mutate_if(is.character,as.numeric)
# Adding a Min and Max rows to the Dataframe to scale the chart
max_min <- data.frame(
N = c(4, 0), NNW = c(4, 0), NW = c(4, 0),
WNW = c(4, 0), W = c(4, 0), WSW = c(4, 0),
SW = c(4, 0), SSW = c(4, 0), S = c(4, 0),SSE = c(4, 0),
SE = c(4, 0), ESE = c(4, 0),E = c(4, 0),
ENE = c(4, 0), NE = c(4, 0), NNE = c(4, 0)
)
rownames(max_min) <- c("Max", "Min")
final_df <- rbind(max_min, final_df)op <- par(mar = c(1, 2, 2, 2))
radarchart(
final_df, axistype = 1,
color = c("#00AFBB","#FC4E07"), pfcol = scales::alpha(c("#00AFBB","#FC4E07"), 0.5),
plwd = 2, plty = 1,cglcol = "grey", cglty = 1, cglwd = 0.8, axislabcol = "grey",
caxislabels = c(0, 1, 2, 3, 4),vlcex = 0.7, vlabels = colnames(final_df),
title = "Average Rainfall in mm by Wind Direction"
)
legend(
x = "bottomleft", legend = rownames(final_df[-c(1,2),]), horiz = TRUE,
bty = "n", pch = 20 , col = c("#00AFBB", "#FC4E07"),
text.col = "black", cex = , pt.cex = 1.5
)
par(op)Criando o modelo
Removendo colunas com muitos dados faltando e depois removendo entradas incompletas.
print(colMeans(is.na(dataset))) # Percentage of missing data in each column## Date Location MinTemp MaxTemp Rainfall
## 0.00000000 0.00000000 0.01020899 0.00866905 0.02241853
## Evaporation Sunshine WindGustDir WindGustSpeed WindDir9am
## 0.43166506 0.48009762 0.07098859 0.07055548 0.07263853
## WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am Humidity3pm
## 0.02906641 0.01214767 0.02105046 0.01824557 0.03098446
## Pressure9am Pressure3pm Cloud9am Cloud3pm Temp9am
## 0.10356799 0.10331363 0.38421559 0.40807095 0.01214767
## Temp3pm RainToday RainTomorrow
## 0.02481094 0.02241853 0.02245978datasetClean <- subset(dataset, select=-c(Evaporation, Sunshine, Cloud9am, Cloud3pm, Date)) # Removing features with high NA precentage
datasetClean <- datasetClean[complete.cases(datasetClean), ] # Removing NAsOne-Hot encoding as features categóricas.
WD3<-dummy(datasetClean$WindDir3pm)
WD9<-dummy(datasetClean$WindDir9am)
WGD<-dummy(datasetClean$WindGustDir)
Loc<-dummy(datasetClean$Location)
colnames(WD3) <- paste("WindDir3pm", colnames(WD3), sep = "_")
colnames(WD9) <- paste("WindDir9am", colnames(WD9), sep = "_")
colnames(WGD) <- paste("WindGustDir", colnames(WGD), sep = "_")
colnames(Loc) <- paste("Location", colnames(Loc), sep = "_")
datasetD <- cbind(datasetClean, WD3)
datasetD <- cbind(datasetD, WD9)
datasetD <- cbind(datasetD, WGD)
datasetD <- cbind(datasetD, Loc)
datasetD$RainToday <- factor(datasetD$RainToday,
levels = c('No', 'Yes'),
labels = c(0, 1))
datasetD$RainTomorrow <- factor(datasetD$RainTomorrow,
levels = c('No', 'Yes'),
labels = c(0, 1))
df<-tibble(datasetD)Separando em datasets de treino e teste.
set.seed(1337)
split = sample.split(df$RainTomorrow, SplitRatio = 0.8)
training_set = subset(df, split == TRUE)
test_set = subset(df, split == FALSE)Treinando o classificador Random Forest.
y <- training_set$RainTomorrow
X <- subset(training_set, select=-c(RainTomorrow))
classifier = randomForest(x = X,
y = y,
ntree = 30)
classifier##
## Call:
## randomForest(x = X, y = y, ntree = 30)
## Type of random forest: classification
## Number of trees: 30
## No. of variables tried at each split: 10
##
## OOB estimate of error rate: 15.53%
## Confusion matrix:
## 0 1 class.error
## 0 66370 3955 0.05623889
## 1 10077 9937 0.50349755Calculando a acurácia do modelo.
# Accuracy Score
y_pred = predict(classifier, newdata = subset(test_set, select=-c(RainTomorrow)))
accuracy = sum(y_pred == test_set$RainTomorrow) / nrow(test_set)
print(accuracy)## [1] 0.8572061Considerando que previsões meteorológicas tem uma acurácia de cerca de 90% para previsões em um período de cinco dias(https://scijinks.gov/forecast-reliability/), a acurácia do modelo de aproximadamente 86% pode ser considerada razoável,
Matriz de confusão.
# Confusion Matrix
cm = table(test_set$RainTomorrow, y_pred)
print(cm)## y_pred
## 0 1
## 0 16774 807
## 1 2418 2586Calculando o precision e recall.
recall = cm[2:2,2]/(cm[2:2,1]+cm[2:2,2])
precision = cm[2:2,2]/(cm[0:1,2]+cm[2:2,2])
cat("Precision:", precision, "\n")## Precision: 0.7621574cat("Recall:", recall)## Recall: 0.5167866Cerca de metade das vezes que o chove modelo prediu que não iria chover. Porém ele acerta a grande maioria das vezes quando prediz que vai ocorrer chuva no dia seguinte.